home *** CD-ROM | disk | FTP | other *** search
/ Graphics Plus / Graphics Plus.iso / general / fractal / kaos.lha / complib / ludcmp.c < prev    next >
Encoding:
C/C++ Source or Header  |  1990-01-27  |  1.3 KB  |  76 lines
  1. /*
  2. ### LU decomposition ###
  3. */
  4.  
  5. #include <math.h>
  6. #define TINY 1.e-20
  7. void ludcmp(a,n,indx,d)
  8. int n,indx[];
  9. double **a,*d;
  10. {
  11.     int i,imax=0,j,k;
  12.     double big,dum,sum,temp;
  13.     double *vv,*dvector(),fabs();
  14.     extern int stop;
  15.  
  16.     vv = dvector(0,n-1);
  17.     *d = 1.0;
  18.     for(i=0;i<n;i++){
  19.         big = 0.0;
  20.         for(j=0;j<n;j++)
  21.             if((temp = fabs(a[i][j])) > big) big = temp;
  22.         if(big == 0.0) {
  23.             printf("ludcmp: Singular matrix\n");
  24.             stop = 1;
  25.             goto done;
  26.         }
  27.         vv[i] = 1.0/big;
  28.     }
  29.     for(j=0;j<n;j++){
  30.         for(i=0;i<j;i++) {
  31.             sum = a[i][j];
  32.             for(k=0;k<i;k++) sum -= a[i][k] * a[k][j];
  33.             a[i][j] = sum;
  34.         }
  35.         big = 0.0;
  36.         for (i=j;i<n;i++) {
  37.             sum = a[i][j];
  38.             for(k=0;k<j;k++) sum -= a[i][k] * a[k][j];
  39.             a[i][j] = sum;
  40.             if((dum = vv[i] * fabs(sum)) >=big) {
  41.                 big = dum;
  42.                 imax = i;
  43.             }
  44.         }
  45.         if(j !=imax) {
  46.             for(k=0;k<n;k++){
  47.                 dum = a[imax][k];
  48.                 a[imax][k] = a[j][k];
  49.                 a[j][k] = dum;
  50.             }
  51.             *d = - (*d);
  52.             vv[imax] = vv[j];
  53.         }
  54.         indx[j] = imax;
  55.         /* modified for our special application */
  56.         /*
  57.         if(a[j][j] == 0) {
  58.             a[j][j] = TINY;
  59.             printf("ludcmp: Singular matrix. Attempt to fix it.\n");
  60.         }
  61.         */
  62.         if(a[j][j] == 0) {
  63.             printf("ludcmp: Singular matrix\n");
  64.             stop = 1;
  65.             goto done;
  66.         }
  67.         if(j != n-1) {
  68.             dum = 1.0 / (a[j][j]);
  69.             for(i=j+1;i<n;i++) a[i][j] *= dum;
  70.         }
  71.     }
  72.  
  73. done:
  74.     (void) free_dvector(vv,0,n-1);
  75. }
  76.